Method and system for determining the position of an object

ABSTRACT

A method and system for determining the location of an object. The system may implement a number of sensors at different locations that are positioned to receive a transmitted or reflected signal from the object. Also disclosed is a system and method of calculating object position based upon the time difference of arrival (TDOA) from each sensor, or the relative time difference of arrival (RTDOA). A known distribution of noise is added to the time of arrival (TOA) prior to calculating the object position.

FIELD OF THE INVENTION

The invention relates to a system and method for determining theposition of an object using a frequency based sensor system (e.g.,radar, sonar, global positioning satellites (GPS), cellular telephony,etc.).

BACKGROUND OF THE INVENTION

Accurately determining the position of moving objects, such as aircraft,missiles, land-vehicles, watercraft, and electronic devices (e.g.,mobile telephones) has been a challenge for sometime. There have beenmany methods developed, including radar, sonar, GPS and other frequencybased techniques.

Typically, object location systems make use of multilaterationtechniques (i.e., the fixing of a position by the time difference ofarrival (TDOA) of a signal at a sensor) to locate the origin of atransmitted signal. TDOA information is typically based on signalarrival-time measurements from an object that transmits, reflects orreceives and re-transmits a signal that is received at receivingsensors. The TDOA measurements and other sensor information is typicallytransmitted to a central processor for processing.

Methods for estimating position using TDOA are known. For example, “ASimple and Efficient Estimator for Hyperbolic Location” by Y. T. Chanand K. C. Ho. (IEEE Transactions on Signal processing, Vol. 42, No. 8,August 1994, pp. 1905–1915), which is incorporated herein by reference,provides techniques for implementing TDOA measurements.

Other methods for determining the location and identification of aobject are also known. For example, U.S. Pat. No. 6,211,811, issued toEvers, discloses a method and apparatus for improving the surveillancecoverage and object identification in a radar based surveillance system.

Another known method is disclosed in U.S. Pat. No. 6,201,499, issued toHawkes et al. As disclosed therein, it is known to provide a method andapparatus for measuring the time difference of arrival between tworeceived signals, each received signal being a time delayed version of atransmitted signal as intercepted at two different sensors where thetransmitted signal is generated by a radio frequency transmitter.

Many drawbacks are present in existing systems and methods. For example,known systems and methods are often inaccurate and costly to implement.Other drawbacks also exist.

SUMMARY OF THE INVENTION

The present invention provides systems and methods for enablingcomputations that can be used to determine the position of a objectusing signal data at multiple sensors. Signal data may comprise any ofTime of Arrival (TOA) data, Time Difference of Arrival (TDOA) data andRelative Time of Arrival (RTOA) data. These computations serve as an aidin the development position locating systems such as next generation airtraffic control radar systems or the like.

Accordingly, there is provided a method for determining the position ofan object in a system comprising multiple sensors arranged at differingheights and including a reference sensor. In some embodiments the methodmay comprise transmitting or reflecting a signal from a object and thendetermining the time of transmission or reflection of the signal fromthe object. The method may then include receiving the transmitted orreflected signal at the multiple sensors and determining the TOA of thesignal at each sensor. A processor-based calculator may calculate aslant range from the object to each sensor and a processor basedcalculator may calculate a position vector of the object.

In accordance with some embodiments of the invention, there is provideda method for determining the position of a object in a system comprisingmultiple sensors arranged at differing heights and including a referencesensor. According to these embodiments, the method may comprisetransmitting or reflecting a signal from a object and determining thetime of transmission or reflection of the signal from the object. Themethod also includes receiving the transmitted or reflected signal atthe multiple sensors and determining the TDOA of the signal at eachsensor. Then, processor based calculators may be used to calculate aslant range from the object to each sensor and to calculate a positionvector for the object.

In accordance with some other embodiments of the invention there isprovided a method for determining the position of an object in a systemcomprising a sensor arranged at a determinable location and a secondarysurveillance device that, at a determinable transmission time, transmitsan interrogator signal that is reflected off the object, or received andretransmitted by the object, to a secondary sensor. In these embodimentsthe method may comprise obtaining a time of arrival for a signalreceived at the sensor and obtaining a secondary time of arrival for thereflected or retransmitted interrogator signal received at the secondarysensor. The method may also include implementing a processor-basedcalculator to calculate a slant range from the object to the sensorbased, at least in part, upon the obtained time of arrival at the sensorand to calculate a secondary slant range from the object to thesecondary sensor based, at least in part, upon the obtained secondarytime of arrival at the secondary sensor. Finally, a processor-basedcalculator may be implemented to determine a position vector based, atleast in part upon the transmission time of the interrogator signal, thecalculated slant range and the calculated secondary slant range.

In accordance with some other embodiments of the invention there isprovided a method for determining the position of an object in a systemcomprising a sensor arranged at a determinable location. In theseembodiments, the method may comprise obtaining a time of arrival for asignal received at the sensor. Then a processor based calculator maycalculate a slant range from the object to the sensor based, at least inpart, upon the obtained time of arrival. Finally a processor basedcalculator may enable a determination of a position vector based, atleast in part, upon the calculated slant range and the location of thesensor.

In accordance with some other embodiments of the invention there isprovided a method for determining the position of an object in a systemcomprising a sensor arranged at a determinable location and a referencesensor. In these embodiments the method may comprise obtaining a timedifference of arrival for a signal received at the sensor with respectto a signal received at the reference sensor and implementing aprocessor based calculator to calculate a slant range from the object tothe sensor based, at least in part, upon the obtained time difference ofarrival. The method may also include enabling a processor basedcalculator to determine a position vector based, at least in part, uponthe calculated slant range and the location of the sensor.

In accordance with some other embodiments of the invention there isprovided a system for determining the position of an object. In theseembodiments the system may comprise a sensor, arranged at a determinablelocation, that obtains a time of arrival for a signal received at thesensor and a secondary surveillance device that, at a determinabletransmission time, transmits an interrogator signal that is reflectedoff the object, or received and retransmitted by the object, to asecondary sensor and obtains a secondary time of arrival for thereflected or retransmitted interrogator signal received at the secondarysensor. The system may also comprise a slant range calculator thatcalculates a slant range from the object to the sensor based, at leastin part, upon the obtained time of arrival at the sensor and a secondaryslant range calculator that calculates a secondary slant range from theobject to the secondary sensor based, at least in part, upon theobtained secondary time of arrival at the secondary sensor. Finally, thesystem may comprise a position vector calculator that determines aposition vector based, at least in part upon the transmission time ofthe interrogator signal, the calculated slant range and the calculatedsecondary slant range.

In accordance with some other embodiments of the invention, there isprovided a system for determining the position of an object. In theseembodiments the system may comprise a sensor arranged at a determinablelocation and a reference sensor, wherein a time difference of arrival isobtained for a signal received at the sensor with respect to a signalreceived at the reference sensor. The system may also comprise a slantrange calculator that calculates a slant range from the object to thesensor based, at least in part, upon the obtained time difference ofarrival and a position vector calculator that determines a positionvector based, at least in part, upon the calculated slant range and thelocation of the sensor.

In accordance with some other embodiments of the invention there isprovided a system for determining the position of an object. In theseembodiments the system may comprise a sensor arranged at a determinablelocation wherein the sensor obtains a time of arrival for a signalreceived at the sensor and a slant range calculator that calculates aslant range from the object to the sensor based, at least in part, uponthe obtained time of arrival. In addition, the system may comprise aposition vector calculator that determines a position vector based, atleast in part, upon the calculated slant range and the location of thesensor.

Other aspects and features of the invention are possible. The followingis a description of some exemplary embodiments of the invention.

BRIEF DESCRIPTION OF THE DRAWINGS

The purpose and advantages of the present invention will be apparent tothose of ordinary skill in the art from the following detaileddescription in conjunction with the appended drawings in which likereference characters are used to indicate like elements.

FIG. 1 is a schematic diagram showing an airborne object and a number ofground-based sensors according to some embodiments of the invention.

FIG. 2 is a three-dimensional schematic diagram showing an object and anumber of sensors according to some embodiments of the invention.

FIG. 3 is a two-dimensional schematic diagram showing an object withfive sensors according to some embodiments of the invention.

FIG. 4 is a plot of TOA computation error contours with three sensorsaccording to one embodiment of the invention.

FIG. 5 is a plot of TOA computation error contours with four sensorsaccording to one embodiment of the invention.

FIG. 6 is a plot of TOA computation error contours with five sensorsaccording to one embodiment of the invention.

FIG. 7 is a three-dimensional TOA schematic diagram showing an objectwith five ground-based sensors according to one embodiment of thepresent invention.

FIG. 8 is a top-view of an error result for a three-dimensionalTOA-based computation with four sensors in accordance with someembodiments of the invention.

FIG. 9 is a side-view of an error result for a three-dimensionalTOA-based computation with four sensors in accordance with someembodiments of the invention.

FIG. 10 is a top-view of an error result for a three-dimensionalTOA-based computation with five sensors in accordance with someembodiments of the invention.

FIG. 11 is a side-view of an error result for a three-dimensionalTOA-based computation with five sensors in accordance with someembodiments of the invention.

FIG. 12 is a top-view of an error result for a three-dimensionalTOA-based computation with six sensors in accordance with someembodiments of the invention.

FIG. 13 is a side-view of an error result for a three-dimensionalTOA-based computation with six sensors in accordance with someembodiments of the invention.

FIG. 14( a) is a plot of TDOA computation error contours with threesensors according to one embodiment of the invention.

FIG. 14( b) is a plot of TDOA computation error contours with foursensors according to one embodiment of the invention.

FIG. 14( c) is a plot of TDOA computation error contours with fivesensors according to one embodiment of the invention.

FIG. 15 is a top-view of an error result for a three-dimensionalCH-TDOA-based computation with four sensors according to someembodiments of the invention.

FIG. 16 is a side-view of an error result for a three-dimensionalCH-TDOA-based computation with four sensors according to someembodiments of the invention.

FIG. 17 is a top-view of an error result for a three-dimensionalCH-TDOA-based computation with five sensors according to someembodiments of the invention.

FIG. 18 is a side-view of an error result for a three-dimensionalCH-TDOA-based computation with five sensors according to someembodiments of the invention.

FIG. 19 is a top-view of an error result for a three-dimensionalCH-TDOA-based computation with six sensors according to some embodimentsof the invention.

FIG. 20 is a side-view of an error result for a three-dimensionalCH-TDOA-based computation with six sensors according to some embodimentsof the invention.

FIG. 21 is a schematic timeline of signals transmitted to a sensor inaccordance with some embodiments of the invention.

FIG. 22 is a schematic flow diagram of a method for computing objectposition adding a known noise distribution to the TOA in accordance withsome embodiments of the invention.

FIG. 23 is a comparison of contour error plots for two-dimensional TOA,CH-TDOA and RTDOA computations with three sensors in accordance withsome embodiments of the invention.

FIG. 24 is a comparison of contour error plots for two-dimensional TOA,CH-TDOA and RTDOA computations with four sensors in accordance with someembodiments of the invention.

FIG. 25 is a comparison of contour error plots for three-dimensionalTOA, CH-TDOA and RTDOA computations with four sensors in accordance withsome embodiments of the invention.

DETAILED DESCRIPTION OF THE INVENTION

The following description is intended to convey a thorough understandingof the invention by providing a number of specific embodiments anddetails. It should be understood, however, that the invention is notlimited to these specific embodiments and details, which are providedfor exemplary purposes only. It should be further understood that onepossessing ordinary skill in the art, in light of known apparatuses andmethods, would appreciate the use of the invention for its intendedpurposes and benefits in any number of alternative embodiments,depending upon specific design and other needs.

The present invention provides a relatively efficient and cost-effectivemethod and apparatus for determining the location of an object. Forexample, FIG. 1 is a schematic representation of an embodiment of theinvention wherein the position vector (i.e., the distance and angularorientation as measured from an origin) of an object 100 may bedetermined via information gathered at a number of sensors 102, 104,106, 108.

As schematically indicated in FIG. 1, object 100 may comprise anaircraft, however, the invention is not so limited. Object 100 maycomprise any object capable of being detected by the appropriate sensors(e.g., 102–108). For example, object 100 may comprise aircraft,missiles, spacecraft (e.g., space shuttles, satellites, etc.), landvehicles, watercraft, or the like. In addition, object 100 may compriseany type of electronic device capable of being detected by appropriatesensors (e.g., 102–108). For example, object 100 may comprise cellulartelephones, paging devices, GPS receivers, wireless radio transmitters,or the like.

Similarly, while four sensors are shown in FIG. 1, the invention is notso limited. Any number of sensors may be used as is appropriate. Inaddition, the type of sensor may depend upon the type of positiondetermining system being employed. For example, position determiningsystems like radar, sonar, GPS, cellular, etc., each employ appropriatetypes of sensors to detect the relevant signals.

It is also possible to practice the claimed invention in otherconfigurations than the one shown in FIG. 1. For example, and asdiscussed in detail below, more or less sensors 102–108 may be employedin a variety of configurations. In addition, it is possible to use amobile object location system wherein the system is contained, at leastin part, in a spacecraft, aircraft, watercraft, or land-vehicle (e.g.,MILSTAR satellite, AWACS or JSTARS aircraft, or similar ship-based orvehicle-based systems).

FIG. 2 is a schematic depiction of an object position detection systemaccording to some embodiments of the invention. As shown in FIG. 2,system 200 may be configured in a 3-dimensional rectangular coordinatesystem. System 200 may comprise a number, M, of sensors arranged atdiffering axis coordinates. As shown in FIG. 2, for some embodiments itis preferable to use six sensors, 205(a), 205(b), 205(c), 205(d),205(e), 205(f), positioned at various locations.

In operation, a transmitted, or reflected signal (not shown) from theobject 203 is received by each of the sensors, 205(a), 205(b), 205(c),205(d), 205(e), 205(f). In some embodiments, the time of arrival (TOA)for the signal from the object to each sensor is recorded for later usein processing. Assuming that the time of transmission, or reflection,from the object is also known, a slant range (i.e., the line of sightdistance between two points, not at the same level relative to aspecific datum) 204(a), 204(b), 204(c), 204(d), 204(e), 204(f) from theobject 203 to each of the sensor 205(a), 205(b), 205(c), 205(d), 205(e),205(f) can be calculated in the following manner:sr _(—)0_(in) =c(t _(—)0_(in) −t _(tgt))  (1)

for in=1, . . . , M where

-   -   sr_(—)0_(in)=Slant range between each sensor in and object,    -   t_(—)0_(in)=Time of arrival,    -   t_(tgt)=Time of transmission or reflection,    -   c=Speed of light.

An example of a slant range can be described with reference to theexample of an airborne radar object (e.g., an airplane flying at highaltitude with respect to a radar antenna). In that example, the slantrange is the hypotenuse of the triangle represented by the altitude ofthe airplane and the distance between the radar antenna and theairplane's ground track.

Because of possible delays and noise contamination of the receivedsignal at each of the sensors 205(a), 205(b), 205(c), 205(d), 205(e),205(f), errors may occur that may skew the accuracy of the TOA of thesignal at each sensor. Therefore, in order to compensate for such anerror, a known noise distribution (e.g., Gaussian noise) may be added tothe TOA, yielding:t _(in) =t _(—)0_(in) +N(σ_(t))  (2)

where N(σ) denotes a Gaussian random number with variance of σ². This isdepicted in FIG. 22, step 2715.

The slant range between sensor in and the object 203, with the addednoise, may be determined by subtracting the TOA of the transmitted orreflected signal from object 203 to the already known time oftransmission or reflection of the signal from the object 203 andmultiplying it with the velocity of the signal (e.g., the speed of light(c)), yielding:r _(in) =c(t _(in) −t _(tgt))  (3)r _(in) =c(t _(—)0_(in) +N(σ_(t))−t _(tgt))r _(in) =c(t _(—)0_(in) −t _(tgt))+cN(σ₁)r _(in) =sr _(—)0_(n) +N(σ_(r)) where σ_(r) =cσ ₁.

A position vector {right arrow over(P)}_(tgt)=[x_(tgt)y_(tgt)z_(tgt)]^(T) of the object 203 maybecalculated based on the TOA of the reflected or transmitted signal ateach sensor t_(in), the slant range from the object to each sensorsr_(in) and the known position of each sensor (i.e., the determinableposition vector {right arrow over (P)}s_(in) of each of the M sensors inthe coordinate system)

${\overset{arrow}{P}s_{in}} = {\begin{bmatrix}x_{in} \\y_{in} \\z_{in}\end{bmatrix}{for}}$in=1, . . . , M.

In some embodiments, a Lorentz inner product is used to calculate objectposition based on TOA information. Following that procedure, let {rightarrow over (v)}₁ and {right arrow over (v)}₂ be a pair of 4-elementvectors:

${{\overset{arrow}{v}}_{1} = \begin{bmatrix}x_{1} \\y_{1} \\z_{1} \\w_{1}\end{bmatrix}},{{\overset{arrow}{v}}_{2} = {\begin{bmatrix}x_{2} \\y_{2} \\z_{2} \\w_{2}\end{bmatrix}.}}$

Then, Lorentz inner product is then defined as

{right arrow over (v)} ₁ ,{right arrow over (v)} ₂

_(L) =x ₁ x ₂ +y ₁ y ₂ +z ₁ z ₂ −w ₁ w ₂.

For each sensor in, the slant range can be expressed the following way:rd _(in) ² =c ²(t _(in) −t _(tgt))² =∥{right arrow over (P)}s _(in)−{right arrow over (P)} _(tgt)∥²=(x _(in) −x _(tgt))²+(y _(in) −y_(tgt))²+(z _(in) −z _(tgt))²,or(x _(in) −x _(tgt))²+(y _(in) −y _(tgt))²+(z _(in) −z _(tgt))² −c ²(t_(in) −t _(tgt))²=0.

Reorganizing the equation gives:2x _(in) x _(tgt)+2y _(in) y _(tgt)+2z _(in) z _(tgt)−2ct _(in) t_(tgt)=(x _(tgt) ² +y _(tgt) ² +z _(tgt) ² −c ² t _(tgt) ²)+(x _(in) ²+y _(in) ² +z _(in) ² −c ² t _(in) ²).

Let λ=x_(tgt) ²+y_(tgt) ²+z_(tgt) ²−c²t_(tgt) ², and organize theequations in=1, . . . , M as a linear system so that2·A·{right arrow over (v)}=λ·I _(—) v+{right arrow over (b)}

where

$\begin{matrix}{{A = \begin{bmatrix}x_{1} & y_{1} & z_{1} & {- {ct}_{1}} \\\vdots & \vdots & \vdots & \vdots \\x_{in} & y_{in} & z_{in} & {- {ct}_{in}} \\\vdots & \vdots & \vdots & \vdots \\x_{M} & y_{M} & z_{M} & {- {ct}_{M}}\end{bmatrix}},} \\{{\overset{arrow}{v} = \begin{bmatrix}x_{tgt} \\y_{tgt} \\z_{tgt} \\{ct}_{tgt}\end{bmatrix}},} \\{{{I\_ v} = \begin{bmatrix}1 \\\vdots \\1 \\\vdots \\1\end{bmatrix}},} \\{\overset{arrow}{b} = \begin{bmatrix}{x_{1}^{2} + y_{1}^{2} + z_{1}^{2} - {c^{2}t_{1}^{2}}} \\\vdots \\{x_{in}^{2} + y_{in}^{2} + z_{in}^{2} - {c^{2}t_{in}^{2}}} \\\vdots \\{x_{M}^{2} + y_{M}^{2} + z_{M}^{2} - {c^{2}t_{M}^{2}}}\end{bmatrix}}\end{matrix}$

Depending on the value of M, the matrix A may or may not be a squarematrix. Assuming that M≧4, the method of least squares is applied bymultiplying both sides with A^(T) so that

$\begin{matrix}{{2A^{T}A\overset{arrow}{v}} = {{\lambda\; A^{T}{I\_ v}} + {A^{T}\overset{arrow}{b}}}} \\{\overset{arrow}{v} = {{\frac{1}{2}( {A^{T}A} )^{- 1}A^{T}{{I\_ v} \cdot \lambda}} + {\frac{1}{2}( {A^{T}A} )^{- 1}A^{T}\overset{arrow}{b}}}} \\{= {{\overset{arrow}{d}\;\lambda} + \overset{arrow}{e}}}\end{matrix}$

where

${\overset{arrow}{d} = {\frac{1}{2}( {A^{T}A} )^{- 1}A^{T}{I\_ v}}},{\overset{arrow}{e} = {\frac{1}{2}( {A^{T}A} )^{- 1}A^{T}{\overset{arrow}{b}.}}}$

Taking the Lorentz inner product of {right arrow over (v)} with itselfgives

$\begin{matrix}{\lambda = {x_{tgt}^{2} + y_{tgt}^{2} + z_{tgt}^{2} - {c^{2}t_{tgt}^{2}}}} \\{= \langle {\overset{arrow}{v},\overset{arrow}{v}} \rangle_{L}} \\{= \langle {{{\overset{arrow}{d} \cdot \lambda} + \overset{arrow}{e}},{{\overset{arrow}{d} \cdot \lambda} + \overset{arrow}{e}}} \rangle_{L}} \\{\lambda = {{\langle {\overset{arrow}{d},\overset{arrow}{d}} \rangle_{L}\lambda^{2}} + {2\langle {\overset{arrow}{d},\overset{arrow}{e}} \rangle_{L}\lambda} + {\langle {\overset{arrow}{e},\overset{arrow}{e}} \rangle_{L}.}}}\end{matrix}$

This results in a quadratic equation with respect to λ:

{right arrow over (d)},{right arrow over (d)}

_(L)λ²+(2

{right arrow over (d)},{right arrow over (e)}

_(L)−1)λ+

{right arrow over (e)},{right arrow over (e)}

_(L)=0.

Solving for λ gives

$\lambda_{\pm} = {\frac{{- ( {{2\langle {\overset{arrow}{d},\overset{arrow}{e}} \rangle_{L}} - 1} )} \pm \sqrt{( {{2\langle {\overset{arrow}{d},\overset{arrow}{e}} \rangle_{L}} - 1} )^{2} - {4\langle {\overset{arrow}{d},\overset{arrow}{d}} \rangle_{L}\langle {\overset{arrow}{e},\overset{arrow}{e}} \rangle_{L}}}}{2\langle {\overset{arrow}{d},\overset{arrow}{d}} \rangle_{L}}.}$

Once λ is calculated, {right arrow over (v)} is obtained by

$\begin{bmatrix}x_{{tgt}_{\pm}} \\y_{{tgt}_{\pm}} \\z_{{tgt}_{\pm}} \\{c \cdot t_{{tgt}_{\pm}}}\end{bmatrix} = {{\overset{arrow}{v}}_{\pm} = {{\overset{arrow}{d} \cdot \lambda_{\pm}} + {\overset{arrow}{e}.}}}$

Now, two possible positions of the object are given by this computation.In order to choose the right solution, the following vectors arecomputed:

${r\_ ch}_{\pm} = {\begin{bmatrix}{{{\overset{arrow}{P}s_{1}} - {\overset{arrow}{P}}_{{tgt}_{\pm}}}} \\\vdots \\{{{\overset{arrow}{P}s_{in}} - {\overset{arrow}{P}}_{{tgt}_{\pm}}}} \\\vdots \\{{{\overset{arrow}{P}s_{M}} - {\overset{arrow}{P}}_{{tgt}_{\pm}}}}\end{bmatrix} - {\begin{bmatrix}{sr}_{1} \\\vdots \\{sr}_{in} \\\vdots \\{sr}_{M}\end{bmatrix}.}}$

The vector with the smallest error norm of r_ch_(±) can be chosen as thecorrect position of the object.

An exemplary calculation of an object position vector is made withreference to FIG. 3. As shown in FIG. 3, a two-dimensional positionlocation system 300, implementing a TOA based calculation, may comprisepositioning four sensors (S₁) 305(b), (S₂) 305(c), (S₃) 305(d), (S₄)305(e), and a reference sensor (S₅) 305(a) in a 20×20 mile 2-dimensionalgrid at locations given by the following sensor position vectors:

$\begin{matrix}{{{\overset{arrow}{P}s_{1}} = \begin{bmatrix}15 \\10\end{bmatrix}},{{\overset{arrow}{P}s_{2}} = \begin{bmatrix}{- 15} \\10\end{bmatrix}},} \\{{{\overset{arrow}{P}s_{3}} = \begin{bmatrix}{- 15} \\{- 10}\end{bmatrix}},{{\overset{arrow}{P}s_{4}} = \begin{bmatrix}15 \\{- 10}\end{bmatrix}},} \\{{\overset{arrow}{P}s_{5}} = {\begin{bmatrix}0 \\0\end{bmatrix}.}}\end{matrix}$

For calculation purposes, each grid node may be considered an objectlocation. For each object position, the position vector calculationcomprises estimating the position of the object for each sensor usingslant range values at each sensor 305(b), 305(c), 305(d), 305(e), tocompute the position of the object.

The above described system and TOA computation may be generalized sothat it can be employed in a three dimensional (3-D) situation. Inaddition, the above system and TOA computation allows for the deploymentof an unlimited number of sensors, so that a more accurate computationof the position of the object may be performed. Then, for each objectposition, the TOA computation may be used to estimate object positionbased on the noise added range data. In addition, a Monte-Carlo approach(or other numerical approximation method) may be employed to compute thestandard deviation of the error between the true position of the objectand the estimated position of the object. The standard deviation ofGaussian error added to the range may be σ_(r)=2.24×10⁻³, or some othersuitable value.

Proceeding with the calculations for each of the sensors sensor 305(a),305(b), 305(c), 305(d), 305(e), the slant range in can be expressed as:sr ² _(in) =c ²(t _(in) −t _(tgt))²  (5)

In terms of the positions of each sensor 305(a), 305(b), 305(c), 305(d),305(e) and the position of the object at each position 303(a), 303(b),303(c), 303(d), slant ranges 307(a), 307(b), 307(c), and 307(d) may bedetermined using the following relationship:sr ² _(in) =∥{right arrow over (P)}s _(in) −{right arrow over (P)}_(tgt)∥²sr ² _(in)=(x _(in) −x _(tgt))²+(y _(in) −y _(tgt))²+(z _(in) −z_(tgt))²  (6)

In other words, equation (6) may be re-written as:

(x_(in)−x_(tgt))²+(y_(in)−y_(tgt))²+(z_(in)−z_(tgt))²−c²(t_(in)−t_(tgt))²=0,and reorganizing the equation yields,2x _(in) x _(tgt)+2y _(in) y _(tgt)+2z _(in) z _(tgt)−2ct _(in) t_(tgt)=(x _(tgt) ² +y _(tgt) ² +z _(tgt) ² −c ² t _(tgt) ²)+(x _(in) ²+y _(in) ² +z _(in) ² −c ² t _(in) ²)  (7)

letting: λ=x_(tgt) ²+y_(tgt) ²+z_(tgt) ²−c²t_(tgt) ² which is a Lorentzinvariant, and organizing equation (7) for in=1, . . . M, above as alinear system, it becomes,2·A·{right arrow over (v)}=λ·I _(—) v+{right arrow over (b)}  (8)

where,

$\begin{matrix}{{A = \begin{bmatrix}x_{1} & y_{1} & z_{1} & {- {ct}_{1}} \\\vdots & \vdots & \vdots & \vdots \\x_{in} & y_{in} & z_{in} & {- {ct}_{in}} \\\vdots & \vdots & \vdots & \vdots \\x_{M} & y_{M} & z_{M} & {- {ct}_{M}}\end{bmatrix}},} \\{{\overset{arrow}{v} = \begin{bmatrix}x_{tgt} \\y_{tgt} \\z_{tgt} \\{ct}_{tgt}\end{bmatrix}},} \\{{{I\_ v} = \begin{bmatrix}1 \\\vdots \\1 \\\vdots \\1\end{bmatrix}},} \\{\overset{arrow}{b} = \begin{bmatrix}{x_{1}^{2} + y_{1}^{2} + z_{1}^{2} - {c^{2}t_{1}^{2}}} \\\vdots \\{x_{in}^{2} + y_{in}^{2} + z_{in}^{2} - {c^{2}t_{in}^{2}}} \\\vdots \\{x_{M}^{2} + y_{M}^{2} + z_{M}^{2} - {c^{2}t_{M}^{2}}}\end{bmatrix}}\end{matrix}$

and depending on the value of M, the matrix A may or may not be a squarematrix. Assuming that M≧4, the method of least squares may be applied bymultiplying the transpose of A(A^(T)) by both sides of equation (8)above, yielding:2A ^(T) A{right arrow over (v)}=λA ^(T) I _(—) v+A ^(T) {right arrowover (b)}  (9)

solving for {right arrow over (v)},{right arrow over (v)}=½(A ^(T) A)⁻¹ A ^(T) I _(—) v·λ+½(A ^(T) A)⁻¹ A^(T) {right arrow over (b)}{right arrow over (v)}={right arrow over (d)}λ+{right arrow over(e)}  (10)where,{right arrow over (d)}=½(A ^(T) A)⁻¹ A ^(T) I _(—) v and {right arrowover (e)}=½(A ^(T) A)⁻¹ A ^(T) {right arrow over (b)}

Taking a Lorentz inner product of {right arrow over (v)} with itselfgives,λ=

{right arrow over (v)},{right arrow over (v)}

_(L)λ=

{right arrow over (d)}·λ+{right arrow over (e)},{right arrow over(d)}·λ+{right arrow over (e)}

_(L)  (11)λ=

{right arrow over (d)},{right arrow over (d)}

_(L)λ²+2

{right arrow over (d)},{right arrow over (e)}

_(L) λ+

{right arrow over (e)},{right arrow over (e)}

_(L), this results in a quadratic equation with respect to λ, yielding

{right arrow over (d)},{right arrow over (d)}

_(L)λ²+(2

{right arrow over (d)},{right arrow over (e)}

_(L)−1)λ+

{right arrow over (e)},{right arrow over (e)}

_(L)=0  (12)

solving for λ yields,

$\begin{matrix}{\lambda_{\pm} = \frac{{- ( {{2\langle {\overset{arrow}{d},\overset{arrow}{e}} \rangle_{L}} - 1} )} \pm {\sqrt{( {{2\langle {\overset{arrow}{d},\overset{arrow}{e}} \rangle_{L}} - 1} )^{2} - {4\langle {\overset{arrow}{d},\overset{arrow}{d}} \rangle_{L}}}\langle {\overset{arrow}{e},\overset{arrow}{e}} \rangle_{L}}}{2\langle {\overset{arrow}{d},\overset{arrow}{d}} \rangle_{L}}} & (13)\end{matrix}$

From the quadratic equation above, λ can be calculated, and hence,{right arrow over (v)} can be obtained by

$\begin{bmatrix}x_{{tgt}_{\pm}} \\y_{{tgt}_{\pm}} \\z_{{tgt}_{\pm}} \\{c \cdot t_{{tgt}_{\pm}}}\end{bmatrix} = {{\overset{arrow}{v}}_{\pm} = {{\overset{arrow}{d} \cdot \lambda_{\pm}} + {\overset{arrow}{e}.}}}$

Because λ yields two values from the solution of equation (13), whichare two possible positions of the object, therefore, in order toascertain the correct position of the object, the following vectors maybe computed:

${r\_ ch}_{\pm} = {\begin{bmatrix}{{{Ps}_{1} - P_{{tgt}_{\pm}}}} \\\vdots \\{{{\overset{arrow}{P}s_{in}} - P_{{tgt}_{\pm}}}} \\\vdots \\{{{\overset{arrow}{P}s_{M}} - P_{{tgt}_{\pm}}}}\end{bmatrix} - \begin{bmatrix}{sr}_{1} \\\vdots \\{sr}_{in} \\\vdots \\{sr}_{M}\end{bmatrix}}$

From the computation above, the vector with the smallest error norm ofr_ch_(±) may be chosen to be the correct position of the object.

FIG. 4 is a contour plot of the standard deviation of the resultingerror for three sensors used in determining the location of an objectaccording to some embodiments of the invention. Likewise, FIG. 5 is acontour plot of the standard deviation of the resulting error for foursensors used in determining the location of an object according to someembodiments of the invention and FIG. 6 is a contour plot of thestandard deviation of the resulting error for five sensors used indetermining the location of an object according to some embodiments ofthe invention.

As can be seen from FIGS. 4–6, the error is relatively small inside thesensor configuration, and it becomes larger for objects positioned awayfrom the center. In particular, FIG. 4 demonstrates that the errorvalues are particularly large along the rays that are collinear to theline segments connecting sensor pairs. One reason for this large errormay be that the discriminant in the quadratic expression of λ is zero ornearly zero, and that may contribute to large errors in the calculation.Another inference that can be surmised from FIGS. 4–6 is that addingmore sensors dramatically improves the accuracy of the positiondetermination.

Another exemplary embodiment of the invention is discussed withreference to FIG. 7. FIG. 7 is a schematic diagram of another embodimentof position location system 400 using a 3-dimensional coordinate systemto implement a TOA computation. As shown, this embodiment may comprisepositioning six sensors 405(a), 405(b), 405(c), 405(d), 405(e), and405(f) in a 20×20 mile grid at the following position vector locations:

$\begin{matrix}{{{\overset{arrow}{P}s_{1}} = \begin{bmatrix}5 \\0 \\0\end{bmatrix}},{{\overset{arrow}{P}s_{2}} = \begin{bmatrix}10 \\10 \\0.5\end{bmatrix}},} \\{{{\overset{arrow}{P}s_{3}} = \begin{bmatrix}{- 8} \\11 \\1\end{bmatrix}},{{\overset{arrow}{P}s_{4}} = \begin{bmatrix}{- 3} \\{- 9} \\1.5\end{bmatrix}},} \\{{{\overset{arrow}{P}s_{5}} = \begin{bmatrix}8 \\{- 8} \\2\end{bmatrix}},{{\overset{arrow}{P}s_{6}} = {\begin{bmatrix}11 \\2 \\2.5\end{bmatrix}.}}}\end{matrix}$discussed above, for calculation purposes, each grid node may beconsidered to be an object position. For each object position, anestimate of the position of the object for each sensor may be performedby using slant range values at each sensor 405(a), 405(b), 405(c),405(d), 405(e), and 405(f). Then, a similar calculation as above is usedto estimate object position based on the noise added range data.

FIGS. 8–13, are schematic representations of results of the calculationfor each grid node presented in terms of a line segment connecting atrue position of the object represented by a circle, and a simulatedposition of the object, represented by an ‘x’ symbol, so that thegreater the segment length connecting the two points, the worse theerror. In other words, where the circle and the x coincide (e.g.,{circle around (x)}) there is relatively small error.

FIG. 8 is a schematic representation of the above described error plotfor a three dimensional case using four sensors and viewed from the top.FIG. 9 is a schematic of the same three dimensional, four sensor plot asFIG. 8, except the view is from the side.

FIG. 10 is a schematic representation of the above described error plotfor a three dimensional case using five sensors and viewed from the top.FIG. 11 is a schematic of the same three dimensional, five sensor plotas FIG. 10, except the view is from the side.

FIG. 12 is a schematic representation of the above described error plotfor a three dimensional case using six sensors and viewed from the top.FIG. 13 is a schematic of the same three dimensional, six sensor plot asFIG. 12, except the view is from the side.

As demonstrated in FIGS. 8–13, the results of the position calculationfor the above described TOA computation gives quite accurate resultsnear the sensor positions. Away from the sensor positions, the resultsof the TOA computation are not as accurate. The trend in values of themean squared error seems to indicate that adding more sensors mayimprove the accuracy, but not by too much because sensor positions alongthe z-axis may be limited by practical considerations.

Some embodiments of the invention may perform object positioncomputation using a TDOA computation. Some methods of performing TDOAcomputations are known. For example, the above noted paper by Chan andHo describes a TDOA calculation.

In contrast to the TOA computation, the TDOA computation uses thedifference of the time of arrival at each sensor with respect to areference sensor to compute the position of a object. Therefore, the setof equations to be solved is somewhat different. Also, the time oftransmission or reflection may not be needed.

First, a description of a Taylor Series TDOA (TS-TDOA) method ofcomputing object position is provided. Initially, the position of theobject is such thatx _(tgt) =x _(c) +Δxy _(tgt) =y _(c) +Δyz _(tgt) =z _(c) +Δz.

Then, the slant range equation for each in can be expanded as a Taylorseries about (x_(c), y_(c), z_(c)) so that

$\begin{matrix}{{sr}_{in} = \sqrt{( {x_{in} - x_{tgt}} )^{2} + ( {y_{in} - y_{tgt}} )^{2} + ( {z_{in} - z_{tgt}} )^{2}}} \\{= \sqrt{( {x_{in} - ( {x_{c} - {\Delta\; x}} )} )^{2} + ( {y_{in} - ( {y_{c} - {\Delta\; y}} )} )^{2} + ( {z_{in} - ( {z_{c} - {\Delta\; z}} )} )^{2}}} \\{\approx {\sqrt{( {x_{in} - x_{c}} )^{2} + ( {y_{in} - y_{c}} )^{2} + ( {z_{in} - z_{c}} )^{2}} +}} \\{\frac{{{- ( {x_{in} - x_{c}} )}\Delta\; x} - {( {y_{in} - y_{c}} )\Delta\; y} - {( {z_{in} - z_{c}} )\Delta\; z}}{\sqrt{( {x_{in} - x_{c}} )^{2} + ( {y_{in} - y_{c}} )^{2} + ( {z_{in} - z_{c}} )^{2}}} + \ldots}\end{matrix}$

Letting sr_(c in)=√{square root over((x_(in)−x_(c))²+(y_(in)−y_(c))²+(z_(in)−z_(c))²)}{square root over((x_(in)−x_(c))²+(y_(in)−y_(c))²+(z_(in)−z_(c))²)}{square root over((x_(in)−x_(c))²+(y_(in)−y_(c))²+(z_(in)−z_(c))²)} and ignoring higherorder terms gives

${sr}_{in} \approx {{sr}_{c\mspace{14mu}{in}} - {\frac{( {x_{in} - x_{c}} )}{{sr}_{c\mspace{14mu}{in}}}\Delta\; x} - {\frac{( {y_{in} - y_{c}} )}{{sr}_{c\mspace{14mu}{in}}}\Delta\; y} - {\frac{( {z_{in} - z_{c}} )}{{sr}_{c\mspace{14mu}{in}}}\Delta\;{z.}}}$

Now, the range difference equation with respect to sensor 1 may bedefined for each in so that

$\begin{matrix}{{sr\_ n}_{{in}{.1}} \equiv {{sr}_{in} - {sr}_{1}}} \\{\approx {{sr}_{c\mspace{14mu}{in}} - {\frac{( {x_{in} - x_{c}} )}{{sr}_{c\mspace{14mu}{in}}}\Delta\; x} - {\frac{( {y_{in} - y_{c}} )}{{sr}_{c\mspace{14mu}{in}}}\Delta\; y} - {\frac{( {z_{in} - z_{c}} )}{{sr}_{c\mspace{14mu}{in}}}\Delta\; z}}} \\{{- {sr}_{c\mspace{14mu} 1}} + {\frac{( {x_{1} - x_{c}} )}{{sr}_{c\mspace{11mu} 1}}\Delta\; x} + {\frac{( {y_{1} - y_{c}} )}{{sr}_{c\mspace{14mu} 1}}\Delta\; y} + {\frac{( {z_{1} - z_{c}} )}{{sr}_{{c\mspace{14mu} 1}\;}}\Delta\; z}} \\{\approx {{sr}_{c\mspace{14mu}{in}} - {sr}_{c\mspace{14mu} 1} + {( {\frac{( {x_{1} - x_{c}} )}{{sr}_{c\mspace{11mu} 1}} - \frac{( {x_{in} - x_{c}} )}{{sr}_{c\mspace{14mu}{in}}}} )\Delta\; x} +}} \\{{( {\frac{( {y_{1} - y_{c}} )}{{sr}_{c\mspace{11mu} 1}} - \frac{( {y_{in} - y_{c}} )}{{sr}_{c\mspace{14mu}{in}}}} )\Delta\; y} + {( {\frac{( {z_{1} - z_{c}} )}{{sr}_{c\mspace{11mu} 1}} - \frac{( {z_{in} - z_{c}} )}{{sr}_{c\mspace{14mu}{in}}}} )\Delta\;{z.}}}\end{matrix}$

Reorganizing the equations in=2, . . . , M as a linear system givesG _(ts) {right arrow over (v)} _(ts) ={right arrow over (h)} _(ts)

where

${G_{ts} = \begin{bmatrix}{\frac{( {x_{1} - x_{c}} )}{{sr}_{c\mspace{14mu} 1}} - \frac{( {x_{2} - x_{c}} )}{{sr}_{c\mspace{14mu} 2}}} & {\frac{( {y_{1} - y_{c}} )}{{sr}_{c\mspace{14mu} 1}} - \frac{( {y_{2} - y_{c}} )}{{sr}_{c\mspace{14mu} 2}}} & {\frac{( {z_{1} - z_{c}} )}{{sr}_{c\mspace{14mu} 1}} - \frac{( {z_{2} - z_{c}} )}{{sr}_{c\mspace{14mu} 2}}} \\\vdots & \vdots & \vdots \\{\frac{( {x_{1} - x_{c}} )}{{sr}_{c\mspace{14mu} 1}} - \frac{( {x_{in} - x_{c}} )}{{sr}_{c\mspace{14mu}{in}}}} & {\frac{( {y_{1} - y_{c}} )}{{sr}_{c\mspace{14mu} 1}} - \frac{( {y_{in} - y_{c}} )}{{sr}_{c\mspace{14mu}{in}}}} & {\frac{( {z_{1} - z_{c}} )}{{sr}_{c\mspace{14mu} 1}} - \frac{( {z_{in} - z_{c}} )}{{sr}_{c\mspace{14mu}{in}}}} \\\vdots & \vdots & \vdots \\{\frac{( {x_{1} - x_{c}} )}{{sr}_{c\mspace{14mu} 1}} - \frac{( {x_{M} - x_{c}} )}{{sr}_{c\mspace{14mu} M}}} & {\frac{( {y_{1} - y_{c}} )}{{sr}_{c\mspace{14mu} 1}} - \frac{( {y_{M} - y_{c}} )}{{sr}_{c\mspace{14mu} M}}} & {\frac{( {z_{1} - z_{c}} )}{{sr}_{c\mspace{14mu} 1}} - \frac{( {z_{M} - z_{c}} )}{{sr}_{c\mspace{14mu} M}}}\end{bmatrix}},{{\overset{\_}{v}}_{ts} = \begin{bmatrix}{\Delta\; x} \\{\Delta\; y} \\{\Delta\; z}\end{bmatrix}},{{\overset{arrow}{h}}_{ts} = {\begin{bmatrix}{{sr\_ n}_{2,1} - ( {{sr}_{c\mspace{14mu} 2} - {sr}_{c\mspace{14mu} 1}} )} \\{{sr\_ n}_{{in},1} - ( {{sr}_{c\mspace{14mu}{in}} - {sr}_{c\mspace{14mu} 1}} )} \\{{sr\_ n}_{M,1} - ( {{sr}_{c\mspace{14mu} M} - {sr}_{c\mspace{14mu} 1}} )}\end{bmatrix}.}}$

With this linear system, an iterative approach may be used to computethe values of (x_(tgt), y_(tgt), z_(tgt)). One of example of aniterative approach is as follows.

The TS TDOA computation begins with range difference values sr_n_(in,1)already known, for in=2, . . . , M. Then, an initial guess on the valuesof the object position x_(c) ⁰, y_(c) ⁰, z_(c) ⁰) may be made with 0 asan index for initial guess of the iteration. Then, the slant range valuesr_(c in) may be calculated for in=1, . . . , M based on the guess ofthe object position. Once sr_(c in) is calculated, the linear systemG_(ts){right arrow over (v)}_(ts)={right arrow over (h)}_(ts) can besolved as a least squares problem for an over-determined system so that{right arrow over (v)} _(ts)=(G _(ts) ^(T) Q ⁻¹ G _(ts))⁻¹ G _(ts) ^(T)Q ⁻¹ {right arrow over (h)} _(ts)

where Q is the covariance matrix of dimension (M−1)×(M−1) and is relatedto the random errors associated with the difference in range values,nominally set to 1 along the diagonal and 0.5 for the off-diagonalelements. Then, a new set of guess values can now be computed from{right arrow over (v)}_(ts)=[Δx Δy Δz]^(T) so that

$\begin{bmatrix}x_{c}^{N + 1} \\y_{c}^{N + 1} \\z_{c}^{N + 1}\end{bmatrix} = {\begin{bmatrix}x_{c}^{N} \\y_{c}^{N} \\z_{c}^{N}\end{bmatrix} + \begin{bmatrix}{\Delta\; x} \\{\Delta y} \\{\Delta\; z}\end{bmatrix}}$

where N is the iteration index. The steps may be repeated with the newguess of the object position by going back and calculating again. If theinitial guess is reasonably accurate, the solution vector shouldconverge to the position of the object (x_(tgt), y_(tgt), z_(tgt)).

One problem with a TS TDOA computation is that it is an iterativemethod. If the initial guess is poor, then the iteration may notconverge to the correct solution. In particular, when computing theobject position using a computer processor the initial guess has to besomewhat close to the correct object position to achieve convergence.This can be inconvenient. In addition, the TS TDOA uses the first orderTaylor series approximation, and there may be cases where theapproximation gives an incorrect solution as well. However, in someembodiments the TS TDOA computation may be implemented as a reference asdiscussed below.

Another TDOA computation that avoids some of the drawbacks of the TSTDOA and can be used to compute object position is based upon the workof Chan and Ho. Herein this computation is called the CH TDOA. First, aset of difference equations in slant range are defined with respect tosensor 1 so that for in=1, . . . , M,sr _(—) n _(in,1) =sr _(in) −sr ₁.Or,sr_(in) =sr _(—) n _(in,1) +sr ₁.

Squaring both sides gives

$\begin{matrix}{{sr}_{in}^{2} = ( {{{sr}_{\_}n_{{in},1}} + {sr}_{1}} )^{2}} \\{= {{sr\_ n}_{{in},1}^{2} + {2{sr\_ n}_{{in},1}{sr}_{1}} + {{sr}_{1}^{2}.}}}\end{matrix}$

On the other hand,

$\begin{matrix}{{sr}_{in}^{2} = {( {x_{in} - x_{tgt}} )^{2} + ( {y_{in} - y_{tgt}} )^{2} + ( {z_{in} - z_{tgt}} )^{2}}} \\{= {( {x_{in}^{2} + y_{in}^{2} + z_{in}^{2}} ) + ( {x_{tgt}^{2} + y_{tgt}^{2} + z_{tgt}^{2}} ) -}} \\{{2x_{in}x_{tgt}} - {2y_{in}y_{tgt}} - {2z_{in}{z_{tgt}.}}}\end{matrix}$

Equating the two gives:sr _(—) n _(in,1) ²+2sr _(—) n _(in,1) sr ₁ +sr ₁ ²=(x _(in) ² +y _(in)² +z _(in) ²)+(x _(tgt) ² +y _(tgt) ² +z _(tgt) ²)−2x _(in) x _(tgt)−2y_(in) y _(tgt)−2z _(in) z _(tgt).

The terms x_(tgt) ², y_(tgt) ², and z_(tgt) ² make the expressionnonlinear. Therefore, a subtraction is applied using the in=1 equationso that for in=2, . . . , M,

sr_n_(in, 1)² + 2sr_n_(in, 1)sr₁ + sr₁² − sr₁² = (x_(in)² + y_(in)² + z_(in)²) + (x_(tgt)² + y_(tgt)² + z_(tgt)²) − 2x_(in)x_(tgt) − 2y_(in)y_(tgt) − 2z_(in)z_(tgt) − ((x₁² + y₁² + z₁²) + (x_(tgt)² + y_(tgt)² + z_(tgt)²) − 2x₁x_(tgt) − 2y₁y_(tgt) − 2z₁z_(tgt))sr_n_(in, 1)² + 2sr_n_(in, 1)sr₁ = −2(x_(in) − x₁)x_(tgt) − 2(y_(in) − y₁)y_(tgt) − 2(z_(in) − z₁)z_(tgt) + (x_(in)² + y_(in)² + z_(in)²) − (x₁² + y₁² + z₁²)

These equations, for in=2, . . . , M, can be recast as a linear systemso thatG _(a) {right arrow over (v)}={right arrow over (h)}

where

${G_{a} = {- \begin{bmatrix}{2( {x_{2} - x_{1}} )} & {2( {y_{2} - y_{1}} )} & {2( {z_{2} - z_{1}} )} & {2{sr\_ n}_{2,1}} \\\vdots & \vdots & \vdots & \vdots \\{2( {x_{in} - x_{1}} )} & {2( {y_{in} - y_{1}} )} & {2( {z_{in} - z_{1}} )} & {2{sr\_ n}_{{in},1}} \\\vdots & \vdots & \vdots & \vdots \\{2( {x_{M} - x_{1}} )} & {2( {y_{M} - y_{1}} )} & {2( {z_{M} - z_{1}} )} & {2{sr\_ n}_{M,1}}\end{bmatrix}}},{\overset{arrow}{v} = \begin{bmatrix}x_{tgt} \\y_{tgt} \\z_{tgt} \\{rd}_{1}\end{bmatrix}},{\overset{arrow}{h} = {\begin{bmatrix}{{sr\_ n}_{2,1}^{2} - ( {x_{2}^{2} + y_{2}^{2} + z_{2}^{2}} ) + ( {x_{1}^{2} + y_{1}^{2} + z_{1}^{2}} )} \\\vdots \\{{sr\_ n}_{{in},1}^{2} - ( {x_{in}^{2} + y_{in}^{2} + z_{in}^{2}} ) + ( {x_{1}^{2} + y_{1}^{2} + z_{1}^{2}} )} \\\vdots \\{{sr\_ n}_{M,1}^{2} - ( {x_{M}^{2} + y_{M}^{2} + z_{M}^{2}} ) + ( {x_{1}^{2} + y_{1}^{2} + z_{1}^{2}} )}\end{bmatrix}.}}$

Assuming that M>3, the linear system can be solved as a least squaresproblem.

First, the error vector is defined as ψ=G_(a){right arrow over(v)}−{right arrow over (h)}. Since Gaussian errors are involved, theleast squares method is applied to minimizeLS=({right arrow over (h)}−G _(a) {right arrow over(v)})^(T)(ψψ^(T))⁻¹({right arrow over (h)}−G _(a) {right arrow over(v)}).

which is solved by{right arrow over (v)}(G _(a) ^(T)Ψ⁻¹ G _(a))⁻¹ G _(a) ^(T)Ψ⁻¹ {rightarrow over (h)}

where

$\begin{matrix}{{{\Psi = {{Exp}( {\psi\;\psi^{T}} )}},}\;} \\{= {{Expectation}\mspace{14mu}{value}\mspace{14mu}{of}\mspace{14mu}\psi\;{\psi^{T}.}}}\end{matrix}$

The covariance matrix Ψ is approximated byΨ=c ² BQB

where

$\begin{matrix}{{Q = \begin{bmatrix}1 & 0.5 & \cdots & \cdots & 0.5 \\0.5 & 1 & 0.5 & \cdots & \vdots \\\vdots & 0.5 & 1 & ⋰ & \vdots \\\vdots & \cdots & ⋰ & ⋰ & 0.5 \\0.5 & \cdots & \cdots & 0.5 & 1\end{bmatrix}},} \\{B = {{{diag}( {{{sr\_}0_{2}},{{sr\_}0_{3}},\ldots,{{sr\_}0_{M}}} )}.}}\end{matrix}$

The matrix Q is the (M−1)×(M−1) covariance matrix oftime-difference-of-arrival, and it is assumed that there is somecorrelation between the time-difference samples as reflected by thenonzero off-diagonal elements. The matrix B is a diagonal matrix withelements (sr_(—)0₂, sr_(—)0₃, . . . , sr_(—)0_(M)). Since B is not knownprior to the computation, it is approximated by solving{right arrow over (v)} ₀=(G _(a) ^(T) Q ⁻¹ G _(a))⁻¹ G _(a) ^(T) Q ⁻¹{right arrow over (h)},

and determining from {right arrow over (v)}₀ the approximate slant rangevalues(sr_(—)0₂, sr_(—)0₃, . . . , sr_(—)0_(M)).

Once Ψ is computed from the approximate value of B and Q, {right arrowover (v)}₁ is computed from{right arrow over (v)} ₁=(G _(a) ^(T)Ψ⁻¹ G _(a))⁻¹ G _(a) ^(T)Ψ⁻¹ {rightarrow over (h)}where{right arrow over (v)}₁=[{tilde over (x)}_(tgt) {tilde over (y)}_(tgt){tilde over (z)}_(tgt) {tilde over (s)}r₁]^(T).

Now, the computed solution {right arrow over (v)}₁ is refined byimposing the condition that{tilde over (s)}r ₁ ² =e ₁ +e ₂ +e ₃.wheree ₁=(x ₁ −x _(tgt))²,e ₂=(y ₁ −y _(tgt))²,e ₃=(z ₁ −z _(tgt))².

This gives a new over-determined system to solve:e ₁=(x ₁ −{tilde over (x)} _(tgt))²,e ₂=(y ₁ −{tilde over (y)} _(tgt))²,e ₃=(z ₁ −{tilde over (z)} _(tgt))²,e ₁ +e ₂ +e ₃ ={tilde over (s)}r ₁ ².

As a linear system,{tilde over (G)} _(a) {tilde over (v)}={tilde over (h)}

where

${{\overset{\sim}{G}}_{a} = \begin{bmatrix}1 & 0 & 0 & 0 \\0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 \\1 & 1 & 1 & 1\end{bmatrix}},{\overset{\sim}{v} = \begin{bmatrix}e_{1} \\e_{2} \\e_{3}\end{bmatrix}},{\overset{arrow}{h} = {\begin{bmatrix}( {x_{1} - {\overset{\sim}{x}}_{tgt}} )^{2} \\( {y_{1} - {\overset{\sim}{y}}_{tgt}} )^{2} \\( {z_{1} - {\overset{\sim}{z}}_{tgt}} )^{2} \\\_^{2} \\{rd}_{1}\end{bmatrix}.}}$

As a least squares problem, the function to be minimized isLS=({tilde over (h)}−{tilde over (G)} _(a) {tilde over (v)})^(T){tildeover (Ψ)}⁻¹({tilde over (h)}−{tilde over (G)} _(a) {tilde over (v)})

where the covariance matrix {tilde over (Ψ)} is approximated by{tilde over (Ψ)}={tilde over (B)}G _(a) ⁻¹ QG _(a) ^(T) ⁻¹ {tilde over(B)}with{tilde over (B)}=diag({tilde over (x)} _(tgt) −x ₁ , {tilde over (y)}_(tgt) −y ₁ , {tilde over (z)} _(tgt) −z ₁ , {tilde over (s)}r ₁).

Thus,

$\begin{matrix}{{\overset{\sim}{v} = ( {{\overset{\sim}{G}}_{a}^{T}{\overset{\sim}{\Psi}}^{- 1}{\overset{\sim}{G}}_{a}} )^{- 1}},{( {{\overset{\sim}{G}}_{a}^{T}{\overset{\sim}{\Psi}}^{- 1}} )\overset{\sim}{h}}} \\{ {{ {= {{\overset{\sim}{( G }}_{a}^{T}( {{\overset{\sim}{B}}^{- 1}G_{a}^{T}Q^{- 1}G_{a}{\overset{\sim}{B}}^{- 1}} ){\overset{\sim}{G}}_{a}}} )^{- 1} \cdot {\overset{\sim}{( G }}_{a}^{T}}( {{\overset{\sim}{B}}^{- 1}G_{a}^{T}Q^{- 1}G_{a}{\overset{\sim}{B}}^{- 1}} )} ){\overset{\sim}{h}.}}\end{matrix}$

Once {tilde over (v)} is calculated, the refined solution is computed by

$\begin{bmatrix}x_{tgt} \\y_{tgt} \\z_{tgt}\end{bmatrix} = {\begin{bmatrix}{x_{1} \pm \sqrt{e_{1}}} \\{y_{1} \pm \sqrt{e_{2}}} \\{z_{1} \pm \sqrt{e_{3}}}\end{bmatrix}.}$

Because of the square root involved, there are up to eight possiblesolutions. In some embodiments, the solution is selected by choosing thesolution closest to the approximate solution ({tilde over (x)}_(tgt),{tilde over (y)}_(tgt), {tilde over (z)}_(tgt)) that was determinedabove.

In a two dimensional (2-D) case with three sensors, the position of theobject can be found without resorting to least squares method. In thiscase, the difference equations can be written

$\begin{matrix}{{\begin{bmatrix}{x_{2} - x_{1}} & {y_{2} - y_{1}} \\{x_{3} - x_{1}} & {y_{3} - y_{1}}\end{bmatrix}\begin{bmatrix}x_{tgt} \\y_{tgt}\end{bmatrix}} = ( {{\begin{bmatrix}{sr\_ n}_{2,1}^{2} \\{sr\_ n}_{3,1}^{2}\end{bmatrix}{sr}_{1}} +} } \\ {\frac{1}{2}\begin{bmatrix}{{sr\_ n}_{2,1}^{2} - ( {x_{2}^{2} + y_{2}^{2}} ) + ( {x_{1}^{2} + y_{1}^{2}} )} \\{{sr\_ n}_{3,1}^{2} - ( {x_{3}^{2} + y_{3}^{2}} ) + ( {x_{1}^{2} + y_{1}^{2}} )}\end{bmatrix}} ) \\{\begin{bmatrix}x_{tgt} \\y_{tgt}\end{bmatrix} = {\begin{bmatrix}{{a_{1} \cdot {sr}_{1}} + b_{1}} \\{{a_{2} \cdot {sr}_{1}} + b_{2}}\end{bmatrix}.}}\end{matrix}$

and substituted intosr ₁ ²=(x _(tgt) −x ₁)²+(y _(tgt) −y ₁)²to giveα·sr ₁ ² +β·sr ₁+γ=0whereα=(a ₁ ² +a ₂ ²−1),β=a ₁(b ₁ −x ₁)+a ₂(b ₂ −y ₁),γ=(b ₁ −x ₁)²+(b ₂ −y ₁)².

Solving for sr₁ gives

${sr}_{1 \pm} = {\frac{{- \beta} \pm \sqrt{\beta^{2} - {4\alpha\;\gamma}}}{2\;\alpha}.}$

From sr₁,

$\begin{bmatrix}x_{{tgt} \pm} \\y_{{tgt} \pm}\end{bmatrix} = {\begin{bmatrix}{{a_{1} \cdot {sr}_{1 \pm}} + b_{1}} \\{{a_{2} \cdot {sr}_{1 \pm}} + b_{2}}\end{bmatrix}.}$

From the two possible solutions, the correct solution can be determinedby calculating the resulting range values and comparing them with theactual data.

In a three dimensional (3-D) case with four sensors, the position of theobject can be found without resorting to least squares method. In thiscase, the difference equations can be written

${{\lbrack \begin{matrix}{x_{2} - x_{1}} & {y_{2} - y_{1}} & {z_{2} - z_{1}} \\{x_{3} - x_{1}} & {y_{3} - y_{1}} & {z_{3} - z_{1}} \\{x_{4} - x_{1}} & {y_{4} - y_{1}} & {z_{4} - z_{1}}\end{matrix} \rbrack\lbrack \begin{matrix}x_{tgt} \\y_{tgt} \\z_{tgt}\end{matrix} \rbrack} = ( {{\lbrack \begin{matrix}{sr\_ n}_{2.1}^{2} \\{sr\_ n}_{3.1}^{2} \\{sr\_ n}_{4.1}^{2}\end{matrix} \rbrack{sr}_{1}} + {\frac{1}{2}\lbrack \begin{matrix}{{sr\_ n}_{2.1}^{2} - ( {x_{2}^{2} + y_{2}^{2} + z_{2}^{2}} ) + ( {x_{1}^{2} + y_{1}^{2} + z_{1}^{2}} )} \\{{sr\_ n}_{3.1}^{2} - ( {x_{3}^{2} + y_{3}^{2} + z_{3}^{2}} ) + ( {x_{1}^{2} + y_{1}^{2} + z_{1}^{2}} )} \\{{sr\_ n}_{4.1}^{2} - ( {x_{4}^{2} + y_{4}^{2} + z_{4}^{2}} ) + ( {x_{1}^{2} + y_{1}^{2} + z_{1}^{2}} )}\end{matrix} \rbrack}} )},{\lbrack \begin{matrix}x_{tgt} \\y_{tgt} \\z_{tgt}\end{matrix} \rbrack = {\lbrack \begin{matrix}{{a_{1} \cdot {sr}_{1}} + b_{1}} \\{{a_{2} \cdot {sr}_{1}} + b_{2}} \\{{a_{3} \cdot {sr}_{1}} + b_{3}}\end{matrix} \rbrack.}}$

and substituted intosr ₁ ²=(x _(tgt) −x ₁)²+(y _(tgt) −y ₁)²+(z _(tgt) −z ₁)²to giveα·sr ₁ ² +β·sr ₁+γ=0whereα=(a ₁ ² +a ₂ ² +a ₃ ²−1),β=a ₁(b ₁ −x ₁)+a ₂(b ₂ −y ₁)+a ₃(b ₃ −z ₁),γ=(b ₁ −x ₁)²+(b ₂ −y ₁)²+(b ₃ −z ₁)².

Solving for sr₁ gives

${sr}_{1 \pm} = {\frac{{- \beta} \pm \sqrt{\beta^{2} - {4\alpha\;\gamma}}}{2\;\alpha}.}$

From sr₁,

$\lbrack \begin{matrix}x_{{tgt} \pm} \\y_{{tgt} \pm} \\z_{{tgt} \pm}\end{matrix} \rbrack = {\lbrack \begin{matrix}{{a_{1} \cdot {sr}_{1 \pm}} + b_{1}} \\{{a_{2} \cdot {sr}_{1 \pm}} + b_{2}} \\{{a_{3} \cdot {sr}_{1 \pm}} + b_{3}}\end{matrix} \rbrack.}$

From the two possible solutions, the correct solution can be determinedby calculating the resulting range values and comparing with the actualdata.

The CH TDOA may be performed for a nominal 2-D case in the same manneras the TOA computation described above. The sensors are positioned atthe same locations as in the 2-D case for the TOA computation.

For this computation, a grid covering a 20 mile×20 mile region andencompassing the sensors is used so that each grid node may beconsidered to be a object position. Then, for each object position, theCH TDOA method described above may be used to estimate object positionbased on the noise added range data. A Monte-Carlo or other numericapproximation approach may be employed to compute the standard deviationof the error between the true position of the object and the estimatedposition of the object. The standard deviation of Gaussian error addedto the range may be σ_(r)=2.24×10⁻³.

In FIGS. 14( a), 14(b) and 14(c), contour plots of the standarddeviation of the resulting errors are shown. As shown in the Figures,the error is small inside the sensor configuration, and it becomeslarger for objects positioned away from the center. In particular, asshown in FIG. 14( a), it is the error values are relatively large alongthe rays that are collinear to the line segments connecting sensorpairs, similar to the results in the TOA computation case for 2-D. Itturns out that along these rays, the discriminant in the quadraticexpression of sr_(1±) above is zero or nearly zero, and that maycontribute to large errors in the CH TDOA calculation.

FIG. 14( b) shows the results for four sensors, there are certainregions inside the sensor configuration where the error increasesdramatically. This result differs greatly from the results of the samesensor configuration with the TOA computation. Examining the matrixG_(a) shows that, if a object is in those regions, either the differenceof the range values or the difference of the components of the sensorpositions are zero or nearly zero. This renders the matrix G_(a) nearlysingular and results in poor estimation of the object position. On theother hand, and as shown in FIG. 14( c), adding an additional sensorremoves these blind spots, and the results appear to be better than thatof the results obtained using the TOA computation.

In a manner similar to the TOA computation, the CH TDOA computation canbe applied to a nominal 3-D case. In the following, the numbers aredimensionless for convenience, but in practice can correspond to anysuitable scale (e.g., miles, kilometers, etc.). The sensors arepositioned at the same locations as those used in the 3-D case for theTOA computation.

As before, a grid covering a 20×20×10 dimensionless region andencompassing the sensors is used so that each grid node is considered tobe a object position. Then, for each object position or grid node, theCH TDOA computation may be used to estimate object position based on thenoise added range data. The standard deviation of Gaussian error addedto the range is σ_(r)=2.24×10⁻³.

FIGS. 15–20 are schematic representations of results of the calculationfor each grid node presented in terms of a line segment connecting atrue position of the object represented by a circle, and a simulatedposition of the object, represented by an ‘x’ symbol, so that thegreater the segment length connecting the two points, the worse theerror. In other words, where the circle and the x coincide (e.g.,{circle around (x)}) there is relatively small error.

FIG. 15 is a schematic representation of the above described error plotfor a three dimensional case using four sensors and viewed from the top.FIG. 16 is a schematic of the same three dimensional, four sensor plotas FIG. 15, except the view is from the side.

FIG. 17 is a schematic representation of the above described error plotfor a three dimensional case using five sensors and viewed from the top.FIG. 18 is a schematic of the same three dimensional, five sensor plotas FIG. 17, except the view is from the side.

FIG. 19 is a schematic representation of the above described error plotfor a three dimensional case using six sensors and viewed from the top.FIG. 20 is a schematic of the same three dimensional, six sensor plot asFIG. 19, except the view is from the side.

In FIGS. 15–20 it is shown that near the sensors, the CH TDOAcomputation gives quite accurate results. Away from the sensors thecomputation is not as accurate.

The plots in FIGS. 17 and 18 indicate that there are certain localizedregions where the G_(a) matrix is once again nearly singular and causingpoor results. Adding one more sensors may dramatically improve theresults as shown in FIGS. 19 and 20.

As with the TOA computation, the trend in the values of mean squarederror for CH TDOA seems to indicate that adding more sensors may improvethe accuracy, but not too much since the sensor positioning along thez-axis is limited.

Comparing the CH TDOA six sensor case in FIGS. 19 and 20 and the TOA sixsensor case in FIGS. 12–13 shows that CH TDOA computation is moreaccurate than the TOA computation especially when multiple sensors areused.

Thus, for some embodiments of the invention, the CH TDOA computationworks better in terms of accuracy when the number of sensors used isfive or more for 2-D and six or more for 3-D. When the number of sensorsis one or two more than the spatial dimension, there may be certainregions or blind spots where the CH TDOA computation fails because theassociated matrix G_(a) becomes nearly singular. In such cases, TOAcomputation appears to work better.

Referring again to FIG. 2, another embodiment of the invention isdescribed wherein an addition method of computing object position isimplemented. As shown in FIG. 2, there are M sensors with differingheights positioned at various locations of the x-y plane. A transmitted,or reflected signal from the object is received by the sensors so thatthe TOA of the signal from the object to each sensor is known. Assumingthat the time of transmission, or reflection, from the object is alsoknown, the slant range from the object to each sensor can be calculated:sr _(—)0_(in) =c(t _(—)0_(in) −t _(tgt))

for in=1, . . . , M where

-   -   sr_(—)0_(in)=Slant range between sensor in and object,    -   t_(—)0_(in)=Time of arrival,    -   t_(tgt)=Time of transmission or reflection,    -   c=Speed of light.

It is likely that the actual times of arrival have errors in them; so, aGaussian noise is added to the time of arrival so that

$\begin{matrix}{{t_{in} = {{{t\_}0_{in}} + {N( \sigma_{t} )}}},} \\{{sr}_{in} = {c( {t_{in} - t_{tgt}} )}} \\{= {c( {{{t\_}0_{in}} + {N( \sigma_{t} )} - t_{tgt}} )}} \\{= {{c( {{{t\_}0_{in}} - t_{tgt}} )} + {{cN}( \sigma_{t} )}}} \\{{{sr}_{in} = {{{sr\_}0_{n}} + {N( \sigma_{r} )}}},}\end{matrix}$

whereσ_(r) =cσ _(t).

The term N(σ) denotes a Gaussian random number with variance of σ².

So, the TDOA computation is to be used to calculate the position vector{right arrow over (P)}_(tgt)=[x_(tgt) y_(tgt) z_(tgt)]^(T) of the objectbased on the following data set:

$\begin{matrix}{{t_{in} = {{Time}\mspace{14mu}{of}\mspace{14mu}{arrival}}},} \\{{{sr}_{in} = {{Slant}\mspace{14mu}{range}\mspace{14mu}{from}\mspace{14mu}{target}\mspace{14mu}{to}\mspace{14mu}{sensor}}},} \\{{{\overset{arrow}{P}s_{in}} = {\lbrack \begin{matrix}x_{in} \\y_{in} \\z_{in}\end{matrix} \rbrack = {{Position}\mspace{14mu}{of}\mspace{14mu}{each}\mspace{14mu}{sensor}}}},}\end{matrix}$

for in=1, . . . , M.

We introduce the following notations for convenience. First, a pair ofnominal vectors is defined as

${\overset{arrow}{x} = \lbrack \begin{matrix}x_{1} \\x_{2} \\x_{3}\end{matrix} \rbrack},{\overset{arrow}{y} = {\lbrack \begin{matrix}y_{1} \\y_{2} \\y_{3}\end{matrix} \rbrack.}}$

Then, an inner product is defined as

{right arrow over (x)}, {right arrow over (y)}

=x ₁ y ₁ +x ₂ y ₂ +x ₃ y ₃.

The norm of a vector is then defined as∥{right arrow over (x)}∥=√{square root over (

{right arrow over (x)}, {right arrow over (y)}

)}.

Using the above notations, for in=2, . . . , M:

$\begin{matrix}{{c( {t_{in} - t_{1}} )} = {c( {t_{in} - t_{tgt} + t_{tgt} - t_{1}} )}} \\{= {{c( {t_{in} - t_{tgt}} )} - {c( {t_{1} - t_{tgt}} )}}} \\{= {{sr}_{in} - {{sr}_{1}.}}}\end{matrix}$

On the other hand,sr _(in) =∥{right arrow over (P)}s _(in) −{right arrow over (P)}_(tgt)∥.

Combining the two expressions givesc(t _(in) −t ₁)=∥{right arrow over (P)}s _(in) ={right arrow over (P)}_(tgt) ∥−∥{right arrow over (P)}s ₁ −{right arrow over (P)} _(tgt)∥.

Multiplying both sides by ∥{right arrow over (P)}s_(in)−{right arrowover (P)}_(tgt)∥+∥{right arrow over (P)}s₁−{right arrow over (P)}_(tgt)∥to remove the square roots gives

$\begin{matrix}{{{c( {t_{in} - t_{1}} )}( {{{{\overset{arrow}{P}s_{in}} - {\overset{arrow}{P}}_{tgt}}} + {{{\overset{arrow}{P}s_{1}} - {\overset{arrow}{P}}_{tgt}}}} )} = {( {{{{\overset{arrow}{P}s_{in}} - {\overset{arrow}{P}}_{tgt}}} - {{{\overset{arrow}{P}s_{1}} - {\overset{arrow}{P}}_{tgt}}}} )( {{{{\overset{arrow}{P}s_{in}} - {\overset{arrow}{P}}_{tgt}}} + {{{\overset{arrow}{P}s_{1}} - {\overset{arrow}{P}}_{tgt}}}} )}} \\{= {{{{\overset{arrow}{P}s_{in}} - {\overset{arrow}{P}}_{tgt}}}^{2} + {{{{\overset{arrow}{P}s_{in}} - {\overset{arrow}{P}}_{tgt}}} \cdot {{{\overset{arrow}{P}s_{1}} - {\overset{arrow}{P}}_{tgt}}}} -}} \\{{{{{\overset{arrow}{P}s_{1}} - {\overset{arrow}{P}}_{tgt}}} \cdot {{{\overset{arrow}{P}s_{in}} - {\overset{arrow}{P}}_{tgt}}}} - {{{\overset{arrow}{P}s_{1}} - {\overset{arrow}{P}}_{tgt}}}^{2}} \\{= {{{{\overset{arrow}{P}s_{in}} - {\overset{arrow}{P}}_{tgt}}}^{2} - {{{\overset{arrow}{P}s_{1}} - {\overset{arrow}{P}}_{tgt}}}^{2}}} \\{= {\langle {{{\overset{arrow}{P}s_{in}} - {\overset{arrow}{P}}_{tgt}},{{\overset{arrow}{P}s_{in}} - {\overset{arrow}{P}}_{tgt}}} \rangle - \langle {{\overset{arrow}{P}s_{1}} - {\overset{arrow}{P}}_{tgt} - {\overset{arrow}{P}s_{1}} - {\overset{arrow}{P}}_{tgt}} \rangle}} \\{= {\langle {{\overset{arrow}{P}s_{in}},{\overset{arrow}{P}s_{in}}} \rangle - {2\langle {{\overset{arrow}{P}s_{1}},{\overset{arrow}{P}}_{tgt}} \rangle} + \langle {{\overset{arrow}{P}s_{tgt}},{\overset{arrow}{P}}_{tgt}} \rangle -}} \\{\langle {{\overset{arrow}{P}s_{1}},{\overset{arrow}{P}s_{1}}} \rangle + {2\langle {{\overset{arrow}{P}s_{1}},{\overset{arrow}{P}}_{tgt}} \rangle} - \langle {{\overset{arrow}{P}s_{tgt}},{\overset{arrow}{P}}_{tgt}} \rangle} \\{= {{{\overset{arrow}{P}s_{in}}}^{2} - {{\overset{arrow}{P}s_{1}}}^{2} + {\langle {{{2\overset{arrow}{P}s_{1}} - {2\overset{arrow}{P}s_{in}}},{\overset{arrow}{P}s_{tgt}}} \rangle.}}}\end{matrix}$

On the other hand,

$\begin{matrix}{{{c( {t_{in} - t_{1}} )}( {{{{\overset{arrow}{P}s_{in}} - {\overset{arrow}{P}}_{tgt}}} + {{{\overset{arrow}{P}s_{1}} - {\overset{arrow}{P}}_{tgt}}}} )} = {{c( {t_{in} - t_{1}} )}( {{c( {t_{in} - t_{tgt}} )} + {c( {t_{1} - t_{tgt}} )}} )}} \\{= {{c^{2}( {t_{in} - t_{1}} )}( {t_{in} + t_{1} - {2t_{tgt}}} )}} \\{= {{{c^{2}( {t_{in} - t_{1}} )}( {t_{in} - t_{1}} )} - {2{c^{2}( {t_{in} - t_{1}} )}{t_{tgt}.}}}}\end{matrix}$

Combining the two gives:∥{right arrow over (P)}s _(in)∥² −∥{right arrow over (P)}s ₁∥²+

2{right arrow over (P)}s ₁ −2{right arrow over (P)}s _(in) ,{right arrowover (P)}s _(tgt)

=c ²(t _(in) −t ₁)(t _(in) +t ₁)−2c ²(t _(in) −t ₁)t _(tgt).

Rearranging it gives:

Equation 1:

2{right arrow over (P)}s ₁−2{right arrow over (P)}s _(in) ,{right arrowover (P)}s _(tgt)

+2c ²(t _(in) −t ₁)t _(tgt) =∥{right arrow over (P)}s ₁∥² −∥{right arrowover (P)}s _(in)∥² +c ²(t _(in) −t ₁)(t _(in) +t ₁).

Now, let the unknown be a vector:

$\overset{arrow}{v} = {\lbrack \begin{matrix}x_{tgt} \\y_{tgt} \\z_{tgt} \\t_{tgt}\end{matrix} \rbrack.}$

Then, Equation 1 can be set up as a linear system for in=2, . . . , M sothat

${{\lbrack \begin{matrix}{{2x_{1}} - {2x_{2}}} & {{2y_{1}} - {2y_{2}}} & {{2z_{1}} - {2z_{2}}} & {2{c( {t_{2} - t_{1}} )}} \\\vdots & \vdots & \vdots & \vdots \\{{2x_{1}} - {2x_{in}}} & {{2y_{1}} - {2y_{in}}} & {{2z_{1}} - {2z_{in}}} & {2{c( {t_{in} - t_{1}} )}} \\\vdots & \vdots & \vdots & \vdots \\{{2x_{1}} - {2x_{M}}} & {{2y_{1}} - {2y_{M}}} & {{2z_{1}} - {2z_{M}}} & {2{c( {t_{M} - t_{1}} )}}\end{matrix} \rbrack\lbrack \begin{matrix}x_{tgt} \\y_{tgt} \\z_{tgt} \\t_{tgt}\end{matrix} \rbrack} = \lbrack \begin{matrix}{{{\overset{arrow}{P}s_{1}}}^{2} - {{\overset{arrow}{P}s_{2}}}^{2} + {{c^{2}( {t_{2} - t_{1}} )}( {t_{2} + t_{1}} )}} \\\vdots \\{{{\overset{arrow}{P}s_{1}}}^{2} - {{\overset{arrow}{P}s_{in}}}^{2} + {{c^{2}( {t_{in} - t_{1}} )}( {t_{in} + t_{1}} )}} \\\vdots \\{{{\overset{arrow}{P}s_{1}}}^{2} - {{\overset{arrow}{P}s_{M}}}^{2} + {{c^{2}( {t_{M} - t_{1}} )}( {t_{M} + t_{1}} )}}\end{matrix} \rbrack},$

orG·{right arrow over (v)}={right arrow over (h)}.

In yet another embodiment of the invention, there is provided a methodfor computing object position using relative time of arrival (RTOA)signals. With reference to FIG. 2, there may be Ns sensors withdiffering heights positioned at various locations of the x-y plane. Inaddition, there may be an additional radar source. For example, asecondary surveillance radar (SSR) may be incorporated as describedbelow.

In FIG. 21, a time history of a signal from the SSR to each of thereceiver sensors is shown. The interrogator signal from the SSR radar istransmitted at t₀ and is reflected by an object at t_(tgt). Thereflected signal is then is received by the sensors so that the time ofarrival of the signal from the object to each sensor is known.

The slant range from the object to each sensor can be calculated as:sr _(in) =c(t _(in) −t _(tgt))

for in=1, . . . , Ns where

-   -   sr_(in)=Slant range between sensor in and object,    -   t_(in)=Time of arrival,    -   t_(tgt)=Time of reflection,    -   c=Speed of light.

The slant range from the object to the SSR is calculated as:sr ₀ =c(t _(tgt) −t ₀).

If the position of the object is represented as {right arrow over(P)}_(tgt)=[x_(tgt) y_(tgt) z_(tgt)]^(T), the slant range between thesensor in and the object issr _(in) =∥{right arrow over (P)}s _(in) −{right arrow over (P)}_(tgt)∥.

So, the proposed RTOA computation is used to calculate the positionvector {right arrow over (P)}_(tgt)=[x_(tgt) y_(tgt) z_(tgt)]^(T) of theobject based on the following data set:

$\begin{matrix}{{t_{0} = {{Time}\mspace{14mu}{of}\mspace{14mu}{transmission}\mspace{14mu}{from}\mspace{14mu}{SSR}}},} \\{{t_{in} = {{Time}\mspace{14mu}{of}\mspace{14mu}{arrival}}},} \\{{{\overset{arrow}{P}s_{in}} = {\lbrack \begin{matrix}x_{in} \\y_{in} \\z_{in}\end{matrix} \rbrack = {{Position}\mspace{14mu}{of}\mspace{14mu}{each}\mspace{14mu}{sensor}}}},} \\{{{Ps}_{SSR} = {\lbrack \begin{matrix}x_{SSR} \\y_{SSR} \\z_{SSR}\end{matrix} \rbrack = {{Position}\mspace{14mu}{of}\mspace{14mu}{SSR}}}},}\end{matrix}$

for in=1, . . . , Ns.

Using the above notations, for in=1, . . . , Ns:

$\begin{matrix}{{c( {t_{in} - t_{0}} )} = {c( {t_{in} - t_{tgt} + t_{tgt} - t_{0}} )}} \\{= {{c( {t_{in} - t_{tgt}} )} + {c( {t_{tgt} - t_{0}} )}}} \\{= {{sr}_{in} + {{sr}_{0}.}}}\end{matrix}$

On the other hand, it is known thatsr _(in) =∥{right arrow over (P)}s _(in) −{right arrow over (P)} _(tgt)∥and sr ₀ =∥{right arrow over (P)}s _(SSR) −{right arrow over (P)}_(tgt)∥.

Combining the two expressions givesc(t _(in) −t ₀)=∥{right arrow over (P)}s _(in) −{right arrow over (P)}_(tgt) ∥+∥{right arrow over (P)}s _(SSR) −{right arrow over (P)}_(tgt)∥,or∥{right arrow over (P)}s _(in) −{right arrow over (P)} _(tgt) ∥=c(t_(in) −t ₀)−∥{right arrow over (P)}s _(SSR) −{right arrow over (P)}_(tgt)∥.

Squaring both sides gives

$\mspace{76mu}{{{{\overset{arrow}{P}s_{in}} - {\overset{arrow}{P}}_{tgt}}}^{2} = ( {{c( {t_{in} - t_{0}} )} - {{{\overset{arrow}{P}s_{SSR}} - {\overset{arrow}{P}}_{tgt}}}} )^{2}}$${{{Ps}_{in}}^{2} - {2\langle {{\overset{arrow}{P}s_{in}},{\overset{arrow}{P}}_{tgt}} \rangle} + {{\overset{arrow}{P}s_{tgt}}}^{2}} = {{c^{2}( {t_{in} - t_{0}} )}^{2} - {2{c( {t_{in} - t_{0}} )}{{{\overset{arrow}{P}s_{SSR}} - {\overset{arrow}{P}}_{tgt}}}} + {{\overset{arrow}{P}s_{SSR}}}^{2} - {2\langle {{\overset{arrow}{P}s_{SSR}},{\overset{arrow}{P}}_{tgt}} \rangle} + {{\overset{arrow}{P}s_{tgt}}}^{2}}$${{{\overset{arrow}{P}s_{in}}}^{2} - {2\langle {{\overset{arrow}{P}s_{in}},{\overset{arrow}{P}}_{tgt}} \rangle}} = {{c^{2}( {t_{in} - t_{0}} )}^{2} - {2{c( {t_{in} - t_{0}} )}{{{\overset{arrow}{P}s_{SSR}} - {\overset{arrow}{P}}_{tgt}}}} + {{\overset{arrow}{P}s_{SSR}}}^{2} - {2{\langle {{\overset{arrow}{P}s_{SSR}},{\overset{arrow}{P}}_{tgt}} \rangle.}}}$

Reorganizing the equation gives

−2{right arrow over (P)}s _(in) +{right arrow over (P)}s _(SSR) ,{rightarrow over (P)}s _(tgt)

+2c(t _(in) −t ₀)∥{right arrow over (P)}s _(SSR) −{right arrow over (P)}_(tgt) ∥=c ²(t _(in) −t ₀)² −∥{right arrow over (P)}s _(in)∥² +∥{rightarrow over (P)}s _(SSR)∥²

Now, define a new variables{tilde over (r)} ₀ =∥{right arrow over (P)}s _(SSR) −{right arrow over(P)} _(tgt)∥.

This then allows the equations for in=1, . . . , Ns to be recast as alinear system of equationsA ₁ {right arrow over (v)} ₁ =b ₁

with the unknown vector defined as:

${\overset{arrow}{v}}_{1} = \lbrack \begin{matrix}{\overset{\sim}{x}}_{tgt} \\{\overset{\sim}{y}}_{tgt} \\{\overset{\sim}{z}}_{tgt} \\{\overset{\sim}{sr}}_{0}\end{matrix} \rbrack$

and

$A_{1} = \lbrack \begin{matrix}{{{- 2}x_{1}} + {2x_{SSR}}} & {{{- 2}y_{1}} + {2y_{SSR}}} & {{{- 2}{z\;}_{1}} + {2\; z_{SSR}}} & {2{c( {t_{1} - t_{0}} )}} \\\vdots & \vdots & \vdots & \vdots \\{{{- 2}x_{in}} + {2x_{SSR}}} & {{{- 2}y_{in}} + {2y_{SSR}}} & {{{- 2}z_{in}} + {2z_{SSR}}} & {2{c( {t_{in} - t_{0}} )}} \\\vdots & \vdots & \vdots & \vdots \\{{{- 2}x_{Ns}} + {2x_{SSR}}} & {{{- 2}y_{Ns}} + {2y_{SSR}}} & {{{- 2}z_{Ns}} + {2z_{SSR}}} & {2{c( {t_{Ns} - t_{0}} )}}\end{matrix} \rbrack$ $b_{1} = {\lbrack \begin{matrix}{{c^{2}( {t_{1} - t_{0}} )}^{2} - {{\overset{arrow}{P}s_{1}}}^{2} + {{\overset{arrow}{P}s_{SSR}}}^{2}} \\\vdots \\{{c^{2}( {t_{in} - t_{0}} )}^{2} - {{\overset{arrow}{P}s_{in}}}^{2} + {{\overset{arrow}{P}s_{SSR}}}^{2}} \\\vdots \\{{c^{2}( {t_{Ns} - t_{0}} )}^{2} - {{\overset{arrow}{P}s_{Ns}}}^{2} + {{\overset{arrow}{P}s_{SSR}}}^{2}}\end{matrix} \rbrack.}$

The method of least squares can be used to solve this linear system sothat{right arrow over (v)} ₁=(A ₁ ^(T) P ⁻¹ A ₁)⁻¹ A ₁ ^(T) P ⁻¹ ·b ₁

where P is the measurement covariance matrix nominally set as identitymatrix of dimension Ns×Ns.

Using Chan and Ho's technique, the covariance matrix H₁ of the errorvector ev₁=A₁{right arrow over (v)}₁−b₁ is approximated from {rightarrow over (v)}₁:H ₁ =c ² C ₁ PC ₁

where

$C_{1} = \lbrack \begin{matrix}{\overset{\sim}{sr}}_{1} & 0 & 0 & 0 & 0 \\0 & ⋰ & 0 & 0 & 0 \\0 & 0 & {\overset{\sim}{sr}}_{in} & 0 & 0 \\0 & 0 & 0 & ⋰ & 0 \\0 & 0 & 0 & \; & {\overset{\sim}{sr}}_{Ns}\end{matrix} \rbrack$${{\overset{\sim}{sr}}_{in} = {{{\overset{arrow}{P}s_{in}} - {\overset{\sim}{P}}_{tgt}}}},{{\overset{\sim}{P}}_{tgt} = {\lbrack \begin{matrix}{\overset{\sim}{x}}_{tgt} \\{\overset{\sim}{y}}_{tgt} \\{\overset{\sim}{z}}_{tgt}\end{matrix} \rbrack.}}$

Now, a more refined solution vector is computed, if desired, by solving{right arrow over (v)} ₁=(A ₁ ^(T) H ₁ ⁻¹ A ₁)⁻¹ A ₁ ^(T) H ₁ ⁻¹ ·b ₁.

Once {right arrow over (v)}₁ is computed, a refinement step similar tothe one employed by Chan and Ho, can be used to obtain a more accuratesolution.

Consider a system:d ₁=(x _(SSR) −{tilde over (x)} _(tgt))²,d ₂=(y _(SSR) −{tilde over (y)} _(tgt))²,d ₃=(z _(SSR) −{tilde over (z)} _(tgt))²,d ₁ +d ₂ +d ₃ ={tilde over (s)}r ₀ ².

As an over-determined linear system,A ₂ {right arrow over (v)} ₂ =b ₂

where

${A_{2} = \lbrack \begin{matrix}1 & 0 & 0 & 0 \\0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 \\1 & 1 & 1 & 1\end{matrix} \rbrack},{{\overset{arrow}{v}}_{1} = \lbrack \begin{matrix}d_{1} \\d_{2} \\d_{3}\end{matrix} \rbrack},{b_{2} = {\lbrack \begin{matrix}( {x_{SSR} - {\overset{\sim}{x}}_{tgt}} )^{2} \\( {y_{SSR} - {\overset{\sim}{y}}_{tgt}} )^{2} \\( {z_{SSR} - {\overset{\sim}{z}}_{tgt}} )^{2} \\{\overset{\sim}{sr}}_{0}^{2}\end{matrix} \rbrack.}}$

As a least squares problem, the function to be minimized isLS=(h ₂ −A ₂ {right arrow over (v)} ₂)^(T) H ₂ ⁻¹(h ₂ −A ₂ {right arrowover (v)} ₂)

where the covariance matrix H₂ of error vector ev₂=A₂{right arrow over(v)}₂−b₂ is approximated byH ₂=4C ₂ F ₁ C ₂ with

$C_{2} = \lbrack \begin{matrix}{x_{SSR} - {\overset{\sim}{x}}_{tgt}} & 0 & 0 & 0 \\0 & {y_{SSR} - {\overset{\sim}{y}}_{tgt}} & 0 & 0 \\0 & 0 & {y_{SSR} - {\overset{\sim}{y}}_{tgt}} & 0 \\0 & 0 & 0 & {\overset{\sim}{sr}}_{0}\end{matrix} \rbrack$

andF ₁=(A ₁ H ₁ ⁻¹ A ₁)⁻¹.

Thus,

$\begin{matrix}{{\overset{arrow}{v}}_{2} = {{( {A_{2}^{T}H_{2}^{- 1}A_{2}} )^{- 1} \cdot ( {A_{2}^{T}H_{2}^{- 1}} )}h_{2}}} \\{= {{( {{A_{2}^{T}( {C_{2}F_{1}C_{2}} )}^{- 1}A_{2}} )^{- 1} \cdot ( {A_{2}^{T}( {C_{2}F_{1}C_{2}} )}^{- 1} )}h_{2}}} \\{= {{( {{A_{2}^{T}( {{C_{2}( {A_{1}H_{1}^{- 1}A_{1}} )}^{- 1}C_{2}} )}^{- 1}A_{2}} )^{- 1} \cdot ( {A_{2}^{T}( {{C_{2}( {A_{1}H_{1}^{- 1}A_{1}} )}^{- 1}C_{2}} )}^{- 1} )}h_{2}}} \\{= {{( {{A_{2}^{T}( {{C_{2}^{- 1}( {A_{1}H_{1}^{- 1}A_{1}} )}C_{2}^{- 1}} )}A_{2}} )^{- 1} \cdot ( {A_{2}^{T}( {{C_{2}^{- 1}( {A_{1}H_{1}^{- 1}A_{1}} )}C_{2}^{- 1}} )} )}h_{2}}}\end{matrix}$

Once {right arrow over (v)}₂ is computed, the refined solution to theobject position is obtained:x _(tgt) −x _(SSR) ±√{square root over (d₁)}y _(tgt) =y _(SSR) ±√{square root over (d₂)}.z _(tgt) =z _(SSR) ±√{square root over (d₃)}

There are eight possible solutions, and the one that is closest to thecoarse solution [{tilde over (x)}_(tgt) {tilde over (y)}_(tgt) {tildeover (z)}_(tgt)]^(T) is selected as the object position.

The functional block diagram of the RTOA computation according to someembodiments is shown in FIG. 22. The steps of the RTOA computation canbe summarized as follows.

After starting the computation at 2710, the position of the receivesensors, the time of arrival of the signal at the sensors, the positionof the SSR radar, and the time of transmission of the interrogatorsignal are acquired as indicated at 2715. As indicated at 2720, a coarseposition of the object is computed by solving{right arrow over (v)} ₁=(A ₁ ^(T) P ⁻¹ A ₁)⁻¹ A ₁ ^(T) P ⁻¹ ·b ₁.

At 2725, the approximated covariance matrix H₁ is computed fromH ₁ =c ² C ₁ PC ₁.

At 2730, the coarse solution is refined by solving{right arrow over (v)} ₁=(A ₁ ^(T) H ₁ ⁻¹ A ₁)⁻¹ A ₁ ^(T) H ₁ ⁻¹ ·b ₁.

At 2735, all the possible refined solutions of the object position arethen computed by solving{right arrow over (v)} ₂=(A ₂ ^(T)(C ₂ ⁻¹(A ₁ H ₁ ⁻¹ A ₁)C ₂ ⁻¹)A₂)⁻¹·(A ₂ ^(T)(C ₂ ⁻¹(A ₁ H ₁ ⁻¹ A ₁)C ₂ ⁻¹))h ₂,

and, at 2740, usingx _(tgt) =x _(SSR) ±√{square root over (e₁)}y _(tgt) =y _(SSR) ±√{square root over (e₂)}.z _(tgt) =z _(SSR) ±√{square root over (e₃)}

At 2745, the refined solutions are compared with the coarse solution,and the one that is closest to the coarse solution is selected as thecomputed object position. At 2750 the computation ends.

Proceeding similarly to the TOA and CH TDOA methods above, the RTOAcomputation my be performed for a 2D case A Monte-Carlo, or othernumerical approximation, approach is used so that the time of arrivaldata and the time of transmit of the interrogation signal have addedGaussian noise error:t _(in) =t _(—)0_(in) +N(σ_(t))

where

-   -   t_(—)0_(in)=Actual time of arrival    -   N(σ_(t))=Random Gaussian error of variance σ_(t) ².

After multiple iterations, the RMS error is computed for eachcomputation and then compared.

The TOA, CH TDOA and RTOA computations are applied to a 2-D case withthree sensors. The sensors are positioned in units of miles at

${{\overset{arrow}{P}s_{1}} = \begin{bmatrix}15 \\10\end{bmatrix}},{{\overset{arrow}{P}s_{2}} = \begin{bmatrix}{- 15} \\10\end{bmatrix}},{{\overset{arrow}{P}s_{3}} = {\begin{bmatrix}{- 15} \\{- 10}\end{bmatrix}.}}$

The SSR is positioned at

${\overset{arrow}{P}s_{SSR}} = {\begin{bmatrix}0 \\0\end{bmatrix}.}$

A grid covering the ±20 mile by ±20 mile region with 101 grid pointsalong each side is used as object positions. For each grid node, thethree computations are applied with 500 iterations to generate a set oferror data so that the RMS error can be computed. The added Gaussianerrors to the time of arrival data have the variance of

$\sigma_{t}^{2} = {( \frac{10\mspace{14mu}{ft}}{c} )^{2}.}$

In FIG. 23, the resulting contours of the RMS errors for the threemultilateration computations are shown. As shown, both the TOA and theTDOA computations have the same RMS errors, which are characterized by aregion of very high errors behind the sensors. It is also clear that theRTOA computation gives much lower RMS error values, which impliesgreater accuracy.

The three multilateration computations are applied to a 2-D case withfour sensors. The sensors are positioned in units of miles at

${{\overset{arrow}{P}s_{1}} = \begin{bmatrix}15 \\10\end{bmatrix}},{{\overset{arrow}{P}s_{2}} = \begin{bmatrix}{- 15} \\10\end{bmatrix}},{{\overset{arrow}{P}s_{3}} = \begin{bmatrix}{- 15} \\{- 10}\end{bmatrix}},{{\overset{arrow}{P}s_{4}} = {\begin{bmatrix}15 \\{- 10}\end{bmatrix}.}}$

The SSR is positioned at

${\overset{arrow}{P}s_{SSR}} = {\begin{bmatrix}0 \\0\end{bmatrix}.}$

The same grid as described above is also used in this case, and theadded Gaussian errors are set to be the same.

In FIG. 24, the resulting contours of the RMS errors for the threemultilateration computations are shown. As shown, the TDOA computationis more accurate compared to the TOA computation in most areas exceptalong the vertical and horizontal in the middle where the RMS errors arevery large. This is likely due to the geometry of the set up in thatalong those two lines the matrix of the linear system is nearlysingular. Similar to the previous case, the RTOA computation gives muchlower RMS error values, which implies greater accuracy.

The three multilateration computations are applied to a 3-D case withfour sensors. The sensors are positioned in units of miles at

${{\overset{arrow}{P}s_{1}} = \begin{bmatrix}15 \\10 \\0\end{bmatrix}},{{\overset{arrow}{P}s_{2}} = \begin{bmatrix}{- 15} \\10 \\0.1\end{bmatrix}},{{\overset{arrow}{P}s_{3}} = \begin{bmatrix}{- 15} \\{- 10} \\0.2\end{bmatrix}},{{\overset{arrow}{P}s_{4}} = {\begin{bmatrix}15 \\{- 10} \\0.3\end{bmatrix}.}}$

The SSR is positioned at

${\overset{arrow}{P}s_{SSR}} = {\begin{bmatrix}0 \\0 \\0\end{bmatrix}.}$

A grid covering the ±20 mile by ±20 mile region with 101 grid pointsalong the x- and y-axis positioned at the z-axis height of 10 miles isused as object positions; and for each grid node, the three computationsare applied with 500 iterations to generate a set of error data so thatthe RMS error can be computed. The added Gaussian errors to the time ofarrival data have the variance of

$\sigma_{t}^{2} = {( \frac{10\mspace{14mu}{ft}}{c} )^{2}.}$

In FIG. 25, the resulting contours of the RMS errors for the threemultilateration computations are shown. As shown, both the TOA and theTDOA computations have the same RMS errors. It is also clear as beforethat the RTOA computation gives much lower RMS error values, whichimplies greater accuracy.

The RTOA computation appears to work better than the previouslydiscussed TOA and TDOA computation in the test cases that are studiedhere. Adding the time of transmission of the interrogator signal fromthe SSR appears to greatly increase the accuracy of the newmultilateration scheme. On the other hand, it is not always clear how toobtain to for some cases. One possible solution is if the return fromthe transponder also includes information about t₀, then the RTOAcomputation could work well.

The invention now being fully described, it will be apparent to one ofordinary skill in the art that many changes and modifications can bemade thereto without departing from the spirit or scope of the inventionas set forth herein. The foregoing describes the preferred embodimentsof the present invention along with a number of possible alternatives.These embodiments, however, are merely for example and the invention isnot restricted thereto. It will be recognized that various materials andmodifications may be employed without departing from the inventiondescribed above, the scope of which is set forth in the followingclaims.

1. A method for determining the position of an object in a systemcomprising a sensor arranged at a determinable location, the methodcomprising: obtaining a time of arrival for a signal received at thesensor calculating a slant range from the object to the sensor based, atleast in part, upon the obtained time of arrival; wherein calculatingthe slant range further comprises: adding a known distribution of noiseto the time of arrival; prior to calculating the slant range; anddetermining a position vector based, at least in part, upon thecalculated slant range and the location of the sensor.
 2. The method ofclaim 1 wherein the time of arrival is obtained from a signaltransmitted from the object.
 3. The method of claim 1 wherein the timeof arrival is obtained from a signal reflected from the object.
 4. Themethod of claim 1 wherein the known distribution of noise comprises aGaussian noise distribution with a variance of σ².
 5. The method ofclaim 1 wherein determining a position vector further comprises:calculating an error norm for each possible position vector solution;and selecting as the object position vector the position vector solutionwith the smallest error norm.